All code is hidden by default, to show code click on code or select ‘show all code’ on the top right of this document.
reqpkg <- c("tidyverse","feather", "DESeq2", "dplyr", "magrittr", "genefilter", "AnnotationHub", "org.Mm.eg.db", "GO.db", "vsn", "pheatmap", "biomaRt", "curl", "RColorBrewer", "viridis", "ggplot2", "ggExtra", "ggpubr", "fgsea", "DT", "ggfortify", "lfda", "rgl", "pca3d")
for (i in reqpkg) {
print(i)
print(packageVersion(i))
suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "tidyverse"
## [1] '1.3.1'
## [1] "feather"
## [1] '0.3.5'
## [1] "DESeq2"
## [1] '1.24.0'
## [1] "dplyr"
## [1] '1.0.6'
## [1] "magrittr"
## [1] '2.0.1'
## [1] "genefilter"
## [1] '1.66.0'
## [1] "AnnotationHub"
## [1] '2.16.1'
## [1] "org.Mm.eg.db"
## [1] '3.8.2'
## [1] "GO.db"
## [1] '3.8.2'
## [1] "vsn"
## [1] '3.52.0'
## [1] "pheatmap"
## [1] '1.0.12'
## [1] "biomaRt"
## [1] '2.40.5'
## [1] "curl"
## [1] '4.3'
## [1] "RColorBrewer"
## [1] '1.1.2'
## [1] "viridis"
## [1] '0.5.1'
## [1] "ggplot2"
## [1] '3.3.3'
## [1] "ggExtra"
## [1] '0.9'
## [1] "ggpubr"
## [1] '0.4.0'
## [1] "fgsea"
## [1] '1.10.1'
## [1] "DT"
## [1] '0.17'
## [1] "ggfortify"
## [1] '0.4.11'
## [1] "lfda"
## [1] '1.1.3'
## [1] "rgl"
## [1] '0.104.16'
## [1] "pca3d"
## [1] '0.10.2'
options(rgl.useNULL = TRUE)
# knitr::knit_hooks$set(webgl = hook_webgl)
theme_set(theme_pubr())
select <- dplyr::select
splits <- list("WTSPF"=1:18, "KOSPF"=19:36, "WTGF"=37:54, "KOGF"=55:72)
group.colors <- c("SPF:WT" = "#103EFB", "SPF:KO" = "#FC2A1C", "GF:WT" ="#128D15", "GF:KO" = "#FD9226")
source("ggbiplot.R")
A = windowsFont("Helvetica LT Pro Light")
ensembl <- useMart("ensembl")
ensemblMouse <- useDataset("mmusculus_gene_ensembl",mart=ensembl)
mouseProteinCodingGenes <- getBM(attributes=c("ensembl_gene_id","external_gene_name","description"), filters='biotype', values=c('protein_coding'), mart=ensemblMouse)
Make DESeq object.
filePath = '../../../'
df <- read.csv(paste0(filePath, 'KF_RNASeq_counts_filtered.txt')) %>%
dplyr::select(-X) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
set_rownames(.$Geneid) %>%
dplyr::select(-Geneid)
gt_list <- rep(c(rep("WT",18), rep("KO",18)),2)
cond_list <- c(rep("SPF",36),rep("GF",36))
time_list <- rep(c(rep("ZT02",3), rep("ZT06",3), rep("ZT10",3), rep("ZT14",3), rep("ZT18",3), rep("ZT22",3)), 4)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list, Time=time_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
dds
## class: DESeqDataSet
## dim: 21847 72
## metadata(1): version
## assays(1): counts
## rownames(21847): ENSMUSG00000051951 ENSMUSG00000025900 ...
## ENSMUSG00000095041 ENSMUSG00000095742
## rowData names(0):
## colnames(72): KF_01 KF_02 ... KF_71 KF_72
## colData names(3): Genotype Condition Time
Normalize data.
vsd <- vst(dds, blind = FALSE)
head(assay(vsd[,1:5]), 3)
## KF_01 KF_02 KF_03 KF_04 KF_05
## ENSMUSG00000051951 4.406910 4.406910 4.406910 4.406910 4.406910
## ENSMUSG00000025900 4.406910 4.406910 4.406910 4.406910 4.406910
## ENSMUSG00000025902 5.946324 5.921394 5.992112 5.503907 5.774604
colData(vsd)
## DataFrame with 72 rows and 4 columns
## Genotype Condition Time sizeFactor
## <factor> <factor> <factor> <numeric>
## KF_01 WT SPF ZT02 1.09295313026657
## KF_02 WT SPF ZT02 1.13261609555455
## KF_03 WT SPF ZT02 1.02495021156392
## KF_04 WT SPF ZT06 1.08800022292013
## KF_05 WT SPF ZT06 1.11995649245623
## ... ... ... ... ...
## KF_68 KO GF ZT18 0.893688439524359
## KF_69 KO GF ZT18 1.02889018870847
## KF_70 KO GF ZT22 0.997131890140167
## KF_71 KO GF ZT22 0.89087483855546
## KF_72 KO GF ZT22 0.970363030607483
# calculate the variance for each gene
rv <- rowVars(assay(vsd))
# select the ntop genes by variance
select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]
# perform a PCA on the data in assay(x) for the selected genes
pca <- prcomp(t(assay(vsd)[select,]))
intgroup <- c("Condition","Genotype")
intgroup.df <- as.data.frame(colData(vsd)[, intgroup, drop=FALSE])
group <- if (length(intgroup) > 1) {
factor(apply( intgroup.df, 1, paste, collapse=":"))
} else {
colData(object)[[intgroup]]
}
percentVar <- pca$sdev^2 / sum(pca$sdev^2)
# assembly the data for the plot
d <- data.frame(PC1=pca$x[,1], PC2=pca$x[,2], PC3=pca$x[,3], group=group, intgroup.df, name=colnames(vsd)) %>%
mutate(colors=c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
plot3d(d[,1:3], pch=16, size=10, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()
plot3d(d[,1:3], pch=16, size=10, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), cex=0.7, pos=3, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()
plot3d(d[,1:3], size=0.1, col = c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), col=c(rep(group.colors[[1]], 18), rep(group.colors[[2]], 18), rep(group.colors[[3]], 18), rep(group.colors[[4]], 18)))
rglwidget()
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=4) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
p1
Variations of the plot.
ggMarginal(p1 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)
ggMarginal(p1 + theme(legend.position = "none"), type="density", groupColour = TRUE)
I also used ggbiplot to plot the PCA. It doesn’t recalculate the PCA, it just plots it differently.
p <- ggbiplot(pca, groups = group, var.axes = F, ellipse = T, circle = T, obs.scale = 0.5) + ggtitle("PCA plot via ggbiplot") + scale_color_manual(values = group.colors)
p$layers[[1]] <- NULL
p + geom_point(aes(color = groups), size=3)
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=4) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
p2
ggMarginal(p2 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)
ggMarginal(p2 + theme(legend.position = "none"), type="density", groupColour = TRUE)
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=4) +
xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
p3
ggMarginal(p3 + theme(legend.position = "none"), type="boxplot", groupFill = TRUE)
ggMarginal(p3 + theme(legend.position = "none"), type="density", groupColour = TRUE)
s = 2
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=s) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=s) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=s) +
xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors)
ggarrange(p1,p2,p3, nrow=1, common.legend = TRUE)
Colors go from lighter to darker.
#DESeq from data set, using data frame for design
dds2 <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Time + Genotype:Condition+Time)
dds2
## class: DESeqDataSet
## dim: 21847 72
## metadata(1): version
## assays(1): counts
## rownames(21847): ENSMUSG00000051951 ENSMUSG00000025900 ...
## ENSMUSG00000095041 ENSMUSG00000095742
## rowData names(0):
## colnames(72): KF_01 KF_02 ... KF_71 KF_72
## colData names(3): Genotype Condition Time
vsd2 <- vst(dds2, blind = FALSE)
head(assay(vsd2[,1:5]), 3)
## KF_01 KF_02 KF_03 KF_04 KF_05
## ENSMUSG00000051951 4.680675 4.680675 4.680675 4.680675 4.680675
## ENSMUSG00000025900 4.680675 4.680675 4.680675 4.680675 4.680675
## ENSMUSG00000025902 6.091255 6.068107 6.133806 5.682343 5.932073
colData(vsd2)
## DataFrame with 72 rows and 4 columns
## Genotype Condition Time sizeFactor
## <factor> <factor> <factor> <numeric>
## KF_01 WT SPF ZT02 1.09295313026657
## KF_02 WT SPF ZT02 1.13261609555455
## KF_03 WT SPF ZT02 1.02495021156392
## KF_04 WT SPF ZT06 1.08800022292013
## KF_05 WT SPF ZT06 1.11995649245623
## ... ... ... ... ...
## KF_68 KO GF ZT18 0.893688439524359
## KF_69 KO GF ZT18 1.02889018870847
## KF_70 KO GF ZT22 0.997131890140167
## KF_71 KO GF ZT22 0.89087483855546
## KF_72 KO GF ZT22 0.970363030607483
rv <- rowVars(assay(vsd2))
select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]
pca2 <- prcomp(t(assay(vsd2)[select,]))
intgroup <- c("Condition","Genotype","Time")
intgroup.df <- as.data.frame(colData(vsd2)[, intgroup, drop=FALSE])
group <- if (length(intgroup) > 1) {
factor(apply( intgroup.df, 1, paste, collapse=":"))
} else {
colData(object)[[intgroup]]
}
percentVar <- pca2$sdev^2 / sum( pca2$sdev^2 )
d <- data.frame(PC1=pca2$x[,1], PC2=pca2$x[,2], PC3=pca2$x[,3], group=group, intgroup.df, name=colnames(vsd2))
group.colors.t <- c("SPF:WT:ZT02" = "#728EFE", "SPF:WT:ZT06" = "#4A6EFC", "SPF:WT:ZT10" = "#224DFC", "SPF:WT:ZT14" = "#0433F1", "SPF:WT:ZT18" = "#032BC9", "SPF:WT:ZT22" = "#0222A1",
"SPF:KO:ZT02" = "#FDA19B", "SPF:KO:ZT06" = "#FD7C72", "SPF:KO:ZT10" = "#FC564A", "SPF:KO:ZT14" = "#FC2A1C", "SPF:KO:ZT18" = "#DD1203", "SPF:KO:ZT22" = "#B50F03",
"GF:WT:ZT02" = "#6FEC71", "GF:WT:ZT06" = "#26E329", "GF:WT:ZT10" = "#17B51A", "GF:WT:ZT14" = "#128D15", "GF:WT:ZT18" = "#0E6C10", "GF:WT:ZT22" = "#0A480B",
"GF:KO:ZT02" = "#FDAD5D", "GF:KO:ZT06" = "#FD9226", "GF:KO:ZT10" = "#FD850D", "GF:KO:ZT14" = "#DE7002", "GF:KO:ZT18" = "#B65C02", "GF:KO:ZT22" = "#8D4701")
group.colors.t2 <- unlist(lapply(group.colors.t, function(x) rep(x, 3)))
group.colors.t
## SPF:WT:ZT02 SPF:WT:ZT06 SPF:WT:ZT10 SPF:WT:ZT14 SPF:WT:ZT18 SPF:WT:ZT22
## "#728EFE" "#4A6EFC" "#224DFC" "#0433F1" "#032BC9" "#0222A1"
## SPF:KO:ZT02 SPF:KO:ZT06 SPF:KO:ZT10 SPF:KO:ZT14 SPF:KO:ZT18 SPF:KO:ZT22
## "#FDA19B" "#FD7C72" "#FC564A" "#FC2A1C" "#DD1203" "#B50F03"
## GF:WT:ZT02 GF:WT:ZT06 GF:WT:ZT10 GF:WT:ZT14 GF:WT:ZT18 GF:WT:ZT22
## "#6FEC71" "#26E329" "#17B51A" "#128D15" "#0E6C10" "#0A480B"
## GF:KO:ZT02 GF:KO:ZT06 GF:KO:ZT10 GF:KO:ZT14 GF:KO:ZT18 GF:KO:ZT22
## "#FDAD5D" "#FD9226" "#FD850D" "#DE7002" "#B65C02" "#8D4701"
show_col(group.colors.t)
plot3d(d[,1:3], pch=16, size=10, col = group.colors.t2)
rglwidget()
plot3d(d[,1:3], pch=16, size=10, col = group.colors.t2)
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), cex=0.7, pos=3, col = group.colors.t2)
rglwidget()
plot3d(d[,1:3], size=0.1, col = group.colors.t2)
text3d(d[,1], d[,2], d[,3], texts=c(rownames(d)), col=group.colors.t2)
rglwidget()
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=4) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none") +
ggtitle("PCA with time gradient")
p1
Variations of the plot.
ggMarginal(p1, type="boxplot", groupFill = TRUE)
ggMarginal(p1, type="density", groupColour = TRUE)
I also used ggbiplot to plot the PCA. It doesn’t recalculate the PCA, it just plots it differently.
p <- ggbiplot(pca, groups = group, var.axes = F, ellipse = T, circle = T, obs.scale = 0.5) + ggtitle("PCA plot via ggbiplot") + scale_color_manual(values = group.colors.t)
p$layers[[1]] <- NULL
p + geom_point(aes(color = groups), size=3)
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=4) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none")
p2
ggMarginal(p2, type="boxplot", groupFill = TRUE)
ggMarginal(p2, type="density", groupColour = TRUE)
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=4) +
xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none")
p3
ggMarginal(p3, type="boxplot", groupFill = TRUE)
ggMarginal(p3, type="density", groupColour = TRUE)
s = 2
p1 <- ggplot(data=d, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=s) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none")
p2 <- ggplot(data=d, aes_string(x="PC1", y="PC3", color="group")) + geom_point(size=s) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none")
p3 <- ggplot(data=d, aes_string(x="PC2", y="PC3", color="group")) + geom_point(size=s) +
xlab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
ylab(paste0("PC3: ",round(percentVar[3] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = group.colors.t) +
theme(legend.position = "none")
ggarrange(p1,p2,p3, nrow=1)
filePath = '../../../data/R/by_time/'
df.PCA <- function(data, intgroup, colors) {
# calculate the variance for each gene
rv <- rowVars(assay(data))
# select the ntop genes by variance
select <- order(rv, decreasing=TRUE)[seq_len(min(500, length(rv)))]
# perform a PCA on the data in assay(x) for the selected genes
pca2 <- prcomp(t(assay(data)[select,]))
intgroup.df <- as.data.frame(colData(data)[, intgroup, drop=FALSE])
group <- if (length(intgroup) > 1) {
factor(apply( intgroup.df, 1, paste, collapse=":"))
} else {
colData(data)[[intgroup]]
}
percentVar <- pca2$sdev^2 / sum( pca2$sdev^2 )
# assembly the data for the plot
d2 <- data.frame(PC1=pca2$x[,1], PC2=pca2$x[,2], group=group, intgroup.df, name=colnames(data))
p <- ggplot(data=d2, aes_string(x="PC1", y="PC2", color="group")) + geom_point(size=3) +
xlab(paste0("PC1: ",round(percentVar[1] * 100),"% variance")) +
ylab(paste0("PC2: ",round(percentVar[2] * 100),"% variance")) +
coord_fixed() +
scale_color_manual(values = colors)
return(p)
}
df <- read_feather(paste0(filePath, 'time_1_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
gt_list <- rep(c(rep("WT",3), rep("KO",3)),2)
cond_list <- c(rep("SPF",6),rep("GF",6))
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p1 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT2")
df <- read_feather(paste0(filePath, 'time_2_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p2 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT6")
df <- read_feather(paste0(filePath, 'time_3_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p3 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT10")
df <- read_feather(paste0(filePath, 'time_4_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p4 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT14")
df <- read_feather(paste0(filePath, 'time_5_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p5 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT18")
df <- read_feather(paste0(filePath, 'time_6_2.feather')) %>%
.[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
as.data.frame() %>%
set_rownames(.$Geneid) %>% dplyr::select(-Geneid)
condition_list <- data.frame(row.names=colnames(df), Genotype=gt_list, Condition=cond_list)
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")
dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
colData = condition_list,
design = ~ Genotype + Condition + Genotype:Condition)
vsd <- vst(dds, blind = FALSE)
p6 <- df.PCA(vsd, c("Condition","Genotype"), group.colors) + ggtitle("ZT22")
ggarrange(p1,p2,p3,p4,p5,p6, common.legend = T)
sessionInfo()
## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19042)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] grid parallel stats4 stats graphics grDevices utils
## [8] datasets methods base
##
## other attached packages:
## [1] scales_1.1.1 plyr_1.8.6
## [3] pca3d_0.10.2 rgl_0.104.16
## [5] lfda_1.1.3 ggfortify_0.4.11
## [7] DT_0.17 fgsea_1.10.1
## [9] Rcpp_1.0.6 ggpubr_0.4.0
## [11] ggExtra_0.9 viridis_0.5.1
## [13] viridisLite_0.3.0 RColorBrewer_1.1-2
## [15] curl_4.3 biomaRt_2.40.5
## [17] pheatmap_1.0.12 vsn_3.52.0
## [19] GO.db_3.8.2 org.Mm.eg.db_3.8.2
## [21] AnnotationDbi_1.46.1 AnnotationHub_2.16.1
## [23] BiocFileCache_1.8.0 dbplyr_2.1.1
## [25] genefilter_1.66.0 magrittr_2.0.1
## [27] DESeq2_1.24.0 SummarizedExperiment_1.14.1
## [29] DelayedArray_0.10.0 BiocParallel_1.18.1
## [31] matrixStats_0.58.0 Biobase_2.44.0
## [33] GenomicRanges_1.36.1 GenomeInfoDb_1.20.0
## [35] IRanges_2.18.3 S4Vectors_0.22.1
## [37] BiocGenerics_0.30.0 feather_0.3.5
## [39] forcats_0.5.1 stringr_1.4.0
## [41] dplyr_1.0.6 purrr_0.3.4
## [43] readr_1.4.0 tidyr_1.1.3
## [45] tibble_3.0.6 ggplot2_3.3.3
## [47] tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] utf8_1.1.4 tidyselect_1.1.0
## [3] RSQLite_2.2.3 htmlwidgets_1.5.3
## [5] munsell_0.5.0 preprocessCore_1.46.0
## [7] miniUI_0.1.1.1 withr_2.4.1
## [9] colorspace_2.0-1 highr_0.8
## [11] knitr_1.31 rstudioapi_0.13
## [13] ggsignif_0.6.0 labeling_0.4.2
## [15] GenomeInfoDbData_1.2.1 bit64_4.0.5
## [17] farver_2.0.3 vctrs_0.3.8
## [19] generics_0.1.0 xfun_0.20
## [21] R6_2.5.0 locfit_1.5-9.4
## [23] manipulateWidget_0.10.1 bitops_1.0-6
## [25] cachem_1.0.1 assertthat_0.2.1
## [27] promises_1.1.1 nnet_7.3-15
## [29] gtable_0.3.0 affy_1.62.0
## [31] klippy_0.0.0.9500 rlang_0.4.10
## [33] splines_3.6.3 rstatix_0.6.0
## [35] broom_0.7.6 checkmate_2.0.0
## [37] BiocManager_1.30.10 yaml_2.2.1
## [39] abind_1.4-5 modelr_0.1.8
## [41] crosstalk_1.1.1 backports_1.2.1
## [43] httpuv_1.5.5 Hmisc_4.4-2
## [45] tools_3.6.3 affyio_1.54.0
## [47] ellipsis_0.3.2 base64enc_0.1-3
## [49] progress_1.2.2 zlibbioc_1.30.0
## [51] RCurl_1.98-1.2 prettyunits_1.1.1
## [53] rpart_4.1-15 cowplot_1.1.1
## [55] haven_2.3.1 cluster_2.1.2
## [57] fs_1.5.0 data.table_1.13.6
## [59] RSpectra_0.16-0 openxlsx_4.2.3
## [61] reprex_2.0.0 hms_1.0.0
## [63] mime_0.9 evaluate_0.14
## [65] xtable_1.8-4 XML_3.99-0.3
## [67] rio_0.5.16 jpeg_0.1-8.1
## [69] readxl_1.3.1 gridExtra_2.3
## [71] compiler_3.6.3 ellipse_0.4.2
## [73] crayon_1.4.1 htmltools_0.5.1.1
## [75] later_1.1.0.1 Formula_1.2-4
## [77] geneplotter_1.62.0 lubridate_1.7.10
## [79] DBI_1.1.1 rappdirs_0.3.2
## [81] Matrix_1.3-2 car_3.0-10
## [83] cli_2.5.0 pkgconfig_2.0.3
## [85] foreign_0.8-75 xml2_1.3.2
## [87] rARPACK_0.11-0 annotate_1.62.0
## [89] webshot_0.5.2 XVector_0.24.0
## [91] rvest_1.0.0 digest_0.6.27
## [93] rmarkdown_2.6 cellranger_1.1.0
## [95] fastmatch_1.1-0 htmlTable_2.1.0
## [97] shiny_1.6.0 lifecycle_1.0.0
## [99] jsonlite_1.7.2 carData_3.0-4
## [101] limma_3.40.6 fansi_0.4.2
## [103] pillar_1.6.1 lattice_0.20-41
## [105] fastmap_1.1.0 httr_1.4.2
## [107] survival_3.2-7 interactiveDisplayBase_1.22.0
## [109] glue_1.4.2 zip_2.1.1
## [111] png_0.1-7 bit_4.0.4
## [113] stringi_1.5.3 blob_1.2.1
## [115] latticeExtra_0.6-29 memoise_2.0.0
smanzoor@uchicago.edu